////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "BaseUtils.h"

ClassImp(BaseUtils);

Int_t BaseUtils::fgFitPrintLevel = 1;
Double_t BaseUtils::fTol = 1;
Int_t BaseUtils::fMaxIter = 100000000;

Int_t BaseUtils::fgFitNStep=0;
Int_t BaseUtils::fgFitNData=0;
Double_t * BaseUtils::fgFitX = 0x0;
Double_t * BaseUtils::fgFitY = 0x0;
Double_t * BaseUtils::fgFitEx = 0x0;
Double_t * BaseUtils::fgFitEy = 0x0;
BaseUtils::FFunc BaseUtils::fgFitFunc = 0x0;

//===================================================================  
//                               Algorithm 
//===================================================================
Bool_t BaseUtils::SolveEq(Double_t &x1, Double_t &x2, const Double_t K1[], const Double_t K2[])
{
  //solve x1, x2 in
  //K10 + K11 x1 + K12 x2 = 0
  //K20 + K21 x1 + K22 x2 = 0

  /*//test
gSystem->Load("../lib/libBaseUtils.so");
double k1[]={7,8,9}
double k2[]={6,3,1}
double x1=-999, x2=-999
BaseUtils::SolveEq(x1,x2,k1,k2)
x1
x2

root [5] x1
(double)(-2.47368421052631593e+00)
root [6] x2
(double)1.42105263157894735e+00
   */

  const Double_t det = K1[1]*K2[2]-K1[2]*K2[1];
  if(TMath::Abs(det)<EPSILON){
    return kFALSE;
  }

  x1 = (K1[2]*K2[0]-K1[0]*K2[2])/det;
  x2 = (K1[0]*K2[1]-K1[1]*K2[0])/det;

  return kTRUE;
}

Double_t BaseUtils::GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2)
{
  //                   ( a0   a1 ) (v0)
  //df^2 = (v0,v1) *   ( a1   a2 ) (v1) = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1
  //f=f(a,b), v0=df/da, v1=df/db
  //f=a/b, v0=f/a, v1=-f/b

  const Double_t df2 = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1;
  if(df2<EPSILON){
    return 0;
  }
  return TMath::Sqrt( df2 );
}

//Numerical Recipes in C, 1997, 3.1 Polynomial Interpolation and Extrapolation, p109
//Neville's algorithm
//error defined here depends on the choice of the route from ya[ns] to y_final
void BaseUtils::Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint)
{
  /*
  for(Int_t ii=0; ii<n; ii++){
    printf("BaseUtils::Polint %10.10e %03d/%03d: %10.10e %10.10e\n", x, ii, n, xa[ii], ya[ii]);
  }
  */

  Double_t c[n];
  Double_t d[n];
  Double_t dift=-999;

  Int_t ns=0;
  Double_t dif = fabs(x-xa[ns]); 
  for (Int_t i=0; i<n;i++) {
    if ( (dift=fabs(x-xa[i])) < dif) { 
      ns=i; 
      dif=dift; 
    }
    c[i]=ya[i]; 
    d[i]=ya[i];

    //to speed up
    /*
    if(kprint) printf("polint =====> x[%d]: %3.10e\ty[%d]: %3.10e\n", i, xa[i], i, ya[i]);
    */
  }

  y = ya[ns];

  //to speed up
  /*
  if(kprint) printf("======> y0: %3.10e at ns=%d\n", y, ns);
  */

  Int_t did=ns-1;
  for (Int_t m=1;m<=n-1;m++){
    for (Int_t i=0; i<n-m; i++){
      Double_t ho=xa[i]-x; 
      Double_t hp=xa[i+m] -x;
      Double_t w = c[i+1]-d[i];
      Double_t den=ho-hp;
      if ( den == 0.0){
        printf("BaseUtils::Polint error (den=ho-hp) == 0.0!! %e %e %e\n", x, xa[i], xa[i+m]); exit(1);
      }

      c[i]=ho*w/den;
      d[i]=hp*w/den; 
      //if(kprint) printf("======> m %d, c[%d]: %5.6e\t d[%d]: %5.6e\n", m,i,c[i],i, d[i]);
    }
    if( m <= n -1-ns ){
      dy = c[ns];

      //to speed up
      /*
      if(kprint) printf("------------> m %d, c[%d]: %5.6e\n", m, ns, dy);
      */
    }
    else {
      dy = d[did];

      //to speed up
      /*
      if(kprint) printf("------------> m %d, d[%d]: %5.6e\n", m, did, dy);
      */

      did--;
    }
    y += dy;

    //to speed up
    /*
    if(kprint) printf("======> y%d: %3.10e\n",m, y);
    */
  }
}

//===================================================================
//                                Functions
//===================================================================
Double_t BaseUtils::BinomialD(const Double_t *xx, const Double_t *par)
{
  return BinomialD(xx[0], par[0], par[1]);
}

Double_t BaseUtils::BinomialD(const Double_t xx, const Double_t nn, const Double_t pp)
{
  if(pp>1 || pp<0){
    return 0;
  }

  const Double_t lng = TMath::LnGamma(nn + 1) - TMath::LnGamma(nn-xx + 1) - TMath::LnGamma(xx + 1);
  const Double_t prob = ProtectedExp(lng)*ProtectedPow(1-pp, nn-xx)*ProtectedPow(pp, xx);

  return prob;
}

Double_t BaseUtils::ProtectedExp(const Double_t xx)
{
  /*
    root [14] exp(709)
    (const double)8.21840746155497238e+307
    root [15] exp(710)
    Error: exp param[0]=710 up:709 low:-inf out of range (tmpfile):1:
  */
  const Double_t maxexparg = 700;
  
  return TMath::Exp(TMath::Min(maxexparg, xx));
}

Double_t BaseUtils::ProtectedPow(const Double_t xx, const Double_t yy)
{
  if(xx<0){
    printf("BaseUtils::ProtectedPow xx<0 %e\n", xx); exit(1);
  }
  //xx^yy = exp(log(xx^yy))=exp(yy log(xx))
  return ProtectedExp(yy*TMath::Log(xx+EPSILON));
}

Double_t BaseUtils::ResolutionGaus(const Double_t *xx, const Double_t *par)
{
  //npar = 3

  //
  //root [0] TMath::Gaus(3, 0, -9,1)
  //(Double_t)(-4.19314697436659056e-02)
  //root [1] TMath::Gaus(3, 0, 9,1)
  //(Double_t)4.19314697436659056e-02
  //
  //!! width has to be abs-ed!!
  return par[0]*TMath::Gaus(xx[0], par[1], par[1]*TMath::Abs(par[2]), kTRUE);
}

Double_t BaseUtils::GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype)
{
  //npar = ntype*npartype

  Double_t sum = 0;

  for(Int_t itype=0; itype<ntype; itype++){
    sum += basefunc(xx, &(par[npartype*itype]));
  }

  return sum;
}

//===================================================================
//                                Fitting 
//===================================================================
TMatrixDSym * BaseUtils::GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[])
{
  // case of fixed params need to take care
  //http://root.cern.ch/root/html/src/TMinuitMinimizer.cxx.html#Ilx.Z

  const Int_t nfree = tmn->GetNumFreePars();
  const Int_t npar = tmn->GetNumPars();

  TMatrixDSym matfree(nfree); 
  tmn->mnemat(matfree.GetMatrixArray(), nfree);

  TMatrixDSym *matall = new TMatrixDSym(npar);

  //r: row, c: column
  Int_t rfree = 0;
  for(Int_t rr=0; rr<npar; rr++){
    if(ifixed[rr]){
      continue;
    }
    Int_t cfree = 0;
    for(Int_t cc=0; cc<npar; cc++){
      if(ifixed[cc]){
        continue;
      }
      (*matall)[rr][cc] = matfree[rfree][cfree];
      cfree++;
    }
    rfree++;
  }

  /*
  //test->
    for(int ii=0; ii<NparBin(); ii++){
    const Int_t id = NparBuffer()+ii;
    printf("ii %d: err %e %e -- %e\n", ii, fLowErrs[id], sqrt((*matall)[id][id]), fLowErrs[id]-sqrt((*matall)[id][id]));
    for(int jj=0; jj<NparBin(); jj++){
    const Int_t i2 = NparBuffer()+jj;
    printf("%e ", (*matall)[id][i2]);
    }
    printf("\n");
    }
    //<-test
  */

  //check consistency!
  const Double_t eps = 1e-1;
  for(Int_t ipar=0; ipar<npar; ipar++){
    const Double_t cov = (*matall)[ipar][ipar];
    const Double_t sigma = TMath::Sqrt(cov);
    if(TMath::Abs(errs[ipar]-sigma)/errs[ipar]>eps){
      printf("BaseUtils::GetCov warning!! --> err-sigma > eps %d: %e %e %e -- %e\n", ipar, errs[ipar], cov, sigma, eps);
    }
  }

  return matall;
}

Int_t BaseUtils::FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //=================================================================================================
  //kml: 
  //0: chi2, 2: binned likelihood
  //=================================================================================================
  
  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }
  //====>

  TMinuit * mnt = new TMinuit(npar);
  switch(kml){
  case 0:
    mnt->SetFCN(MinFCNChisquare);
    break;
  case 2:
    mnt->SetFCN(MinFCNBinnedLikelihood);
    break;
  default:
    printf("BaseUtils::FitKernel bad kml %d\n", kml);exit(1);
  }
  mnt->SetPrintLevel(fgFitPrintLevel);
    
  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //SET ERRordef<up>
  //Sets the value of UP (default value= 1.), defining parameter errors. Minuit defines parameter errors as the change in parameter value required to change the function value by UP. Normally, for chisquared fits UP=1, and for negative log likelihood, UP=0.5. 
  const Double_t errdef = kml? 0.5 : 1;
  mnt->SetErrorDef(errdef);
    
  Int_t nfree = 0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    //1e-3*par seems too small
    const Double_t ss= TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*6e-1 : 1e-3;//important if limits are set
    //1e-2; //restore to old values for TRD performance /u/xlu/.task/dedx/drawSummary in /u/xlu/.task/dedx/AllDataSummary/hist/draw_UpdatedStyle_Comb
    //
    const Double_t ll=(lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hh=(highs && highs[ipar]>-1e9)?highs[ipar]:0;
    
    mnt->DefineParameter(ipar, Form("p%d",ipar), start, ss, ll,hh);
    nfree++;
    if(pfix && pfix[ipar]){
      mnt->FixParameter(ipar);
      nfree--;
    }
  }

  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //MIGrad[maxcalls] [tolerance]
  //Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

  Double_t arg[2]={fMaxIter, fTol};
  Int_t migradflag = -999;
  mnt->mnexcm("MIGRAD", arg, 2, migradflag);

  //get par right after MIGRAD before HESSE modifies them, initial pars not changed
  Double_t migp[10];
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<npar;ipar++){
    mnt->GetParameter(ipar,migp[ipar], dummy);
  }

  Int_t hesseflag = -999;
  mnt->mnexcm("HESSE", arg ,1, hesseflag);

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("BaseUtils::FitKernel fitting error !! %d %d keep raw par and set err to -999!!!\n", migradflag, hesseflag);
    if(errs){
      for(Int_t ipar=0;ipar<npar;ipar++){
        errs[ipar]=-999;
      }
    }
    if(chi){
      chi[0]=chi[1]=-999;
    }
  }
  else{
    for(Int_t ipar=0;ipar<npar;ipar++){
      if(errs) {
        mnt->GetParameter(ipar, dummy,errs[ipar]);
      }

      pars[ipar]=migp[ipar];
    }
    if(chi){
      chi[0] = GetFCNChisquare(pars);
      if(nfree!=mnt->GetNumFreePars()){
        printf("BaseUtils::FitKernel nfree!=mnt->GetNumFreePars() %d %d\n", nfree, mnt->GetNumFreePars());exit(1);
      }
      chi[1] = -nfree;//need to add nx outside FitKernel
    }
  }

  if(cov){
    mnt->mnemat(cov,npar);
  }

  delete mnt;

  //<===
 
  if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 
 
  //test->
  /*
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
  */
  //<-test

  return kfail;
}

Int_t BaseUtils::BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //=================================================================================================
  //binned likelihood fit, kml=1
  //=================================================================================================
  
  //keep empty bins
  const Double_t thres = -1;
  return BinnedFitHist(kTRUE, hh, ffunc, npar, pars, errs, 0x0, thres, kFALSE, pfix, cov, lows, highs);
}

Int_t BaseUtils::BinnedFitHist(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //=================================================================================================
  //fitting hist, either chisquare (kml=kFALSE) or likelihood (kml=kTRUE)
  //=================================================================================================

  const Int_t nmax = hh->GetNbinsX();
  Double_t xdata[nmax];
  Double_t ydata[nmax];
  Double_t xerr[nmax];
  Double_t yerr[nmax];

  Int_t nx=0;
  for(Int_t ii=hh->GetXaxis()->GetFirst(); ii<=hh->GetXaxis()->GetLast(); ii++){
    const Double_t ey = hh->GetBinError(ii);
    const Double_t yy = hh->GetBinContent(ii);
    //keep all for kml==kTRUE
    //very important for BinnedLikelihood!! only keeping empty bins can the result be good
    if(!kml){
      //skip empty bin
      if(ey<EPSILON){
        if(yy>EPSILON){
          printf("BaseUtilsHist::ChisquareFit error! %d %e %e\n", ii, ey, yy); exit(1);
        }
        continue;
      }

      //skip low statstics bin
      if(yy<thres){
        continue;
      }
    }
    const Double_t xx = hh->GetBinCenter(ii);
    const Double_t ex = hh->GetBinWidth(ii)/TMath::Sqrt(12);

    xdata[nx]=xx;
    ydata[nx]=yy;
    xerr[nx]=ex;
    yerr[nx]=ey;
    nx++;
  }

  return BinnedFitLevel0(kml, nx, xdata, ydata, kxerr?xerr:0x0, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t BaseUtils::ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //=================================================================================================
  //chisquare fit with kml = 0
  //=================================================================================================
  return BinnedFitLevel0(kFALSE, nx, xdata, ydata, xerr, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t BaseUtils::BinnedFitLevel0(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //=================================================================================================
  //routine for binned fit, either maximum likelihood or chisquare fit, from arrays X and Y
  //=================================================================================================

  if(nx<=0){
    printf("BaseUtils::BinnedFitLevel0 error nx<=0 %d\n", nx); exit(1);
  }

  fgFitNStep=0;
  fgFitNData = nx;
  fgFitX = new Double_t[nx];
  fgFitY = new Double_t[nx];
  fgFitEx = new Double_t[nx];
  fgFitEy = new Double_t[nx];
  for(Int_t ix=0; ix<fgFitNData; ix++){
    fgFitX[ix]=xdata[ix];
    fgFitY[ix]=ydata[ix];
    fgFitEx[ix]=xerr?xerr[ix]:0;
    fgFitEy[ix]=yerr[ix];
    if(fgFitPrintLevel>=0){
      printf("BaseUtils::BinnedFitLevel0 %d: %e %e %e %e\n", ix, fgFitX[ix], fgFitY[ix], fgFitEx[ix], fgFitEy[ix]);
    }
  }
  fgFitFunc = ffunc;

  const Int_t kfail = FitKernel(kml?2:0, npar, pars, errs, chi, pfix, cov, lows, highs);

  if(chi){
    chi[1] += nx;
    printf("\tBaseUtils::BinnedFitLevel0 Chi2/NDOF %.1f/%.0f = %.1f\n\n", chi[0], chi[1], chi[1]? chi[0]/chi[1] : -999);
  }

  delete fgFitX;
  delete fgFitY;
  delete fgFitEx;
  delete fgFitEy;

  return kfail;
}

Double_t BaseUtils::GetFCNChisquare(const Double_t pars[])
{
  //=================================================================================================
  //chi2
  //=================================================================================================
  
  //const Int_t ktest = 1000000;//207;//4165;

  Double_t f = 0;

  /*
  if(fgFitNStep>=ktest){
    for(Int_t ii=0; ii<npar; ii++){
      printf("in MinFCN pars %d %e\n", ii, pars[ii]);
    }
  }
  */

  for(Int_t idata=0; idata<fgFitNData; idata++){
    //if(fgFitNStep>=ktest) printf("test %d %e %e %e %e %e %e\n", idata, fgFitX[idata], pars[0], pars[1], pars[2], pars[3], pars[4]);

    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t dy = fgFitY[idata] - evaly;
    
    Double_t dfdx=0;
    //calculated df/dx if necessary
    if( fgFitEx[idata]!=0){
      const Double_t derix0=fgFitX[idata]-fgFitEx[idata];
      const Double_t derix1=fgFitX[idata]+fgFitEx[idata];
      const Double_t deriy0=fgFitFunc(&derix0, pars);
      const Double_t deriy1=fgFitFunc(&derix1, pars);
      dfdx = (deriy1-deriy0)/(derix1-derix0);
    }
    f += dy*dy/(fgFitEy[idata]*fgFitEy[idata] + dfdx*fgFitEx[idata]*dfdx*fgFitEx[idata]);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e data %e f %e dy %e ey %e\n", idata, pars[0], pars[1], fgFitX[idata], f, dy, fgFitEy[idata]);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  fgFitNStep++;

  return f;
}

void BaseUtils::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //=================================================================================================
  //function to be minimized in chi2 fit
  //=================================================================================================
  f = GetFCNChisquare(pars);
}

Double_t BaseUtils::GetFCNBinnedLikelihood(const Double_t pars[])
{
  //=================================================================================================
  //binned likelihood
  //=================================================================================================

  //const Int_t ktest = 0;//1000000;//207;//4165;

  //the normalization is important for the error estimation, but not the mean!
  //the norm should be the same for all bins, so just get it from first non-0 bin
  //no problem with bin-width, perfect!
  Double_t norm = 0;
  for(Int_t idata=0; idata<fgFitNData; idata++){
    if(fgFitEy[idata]>EPSILON){
      if(fgFitY[idata]<EPSILON){
        printf("BaseUtils::GetFCNBinnedLikelihood error y<eps, while ey>eps!! %e %e %e", EPSILON,fgFitEy[idata],fgFitY[idata] );exit(1);
      }
      //y = kN, ey = k sqrt(N), y/ey = sqrt(N)
      //k include absolute count and bin width
      const Double_t trueCount = TMath::Power(fgFitY[idata]/fgFitEy[idata],2);
      norm = trueCount/fgFitY[idata];
      break;
    }
  }
  if(norm<EPSILON){
    printf("BaseUtils::GetFCNBinnedLikelihood error norm 0 !! %e\n", norm); exit(1);
  }

  Double_t f = 0;

  for(Int_t idata=0; idata<fgFitNData; idata++){
    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t prob = TMath::Poisson(fgFitY[idata]*norm, evaly*norm);
    f += -TMath::Log( prob>EPSILON ? prob : EPSILON);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e p2 %e x %e y %e y0 %e prob %e f %e\n", idata, pars[0], pars[1], pars[2], fgFitX[idata], fgFitY[idata], evaly, prob, f);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  //if(fgFitNStep>=ktest) printf("test\n"); exit(1);

  fgFitNStep++;
  return f;
}

void BaseUtils::MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //=================================================================================================
  //function to be minimized in binned likelihood fit
  //=================================================================================================
  f = GetFCNBinnedLikelihood(pars);
}

//===================================================================
//                           Histograms
//===================================================================
Double_t BaseUtils::FindYield(const TH2D *hh, const Double_t target)
{
  /*
gSystem->Load("../lib/libBaseUtils.so")
hh=horiginal
new TCanvas; hh->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()
hhy=BaseUtils::GetHyield(hh, "xxx");
new TCanvas; hhy->Draw(); gStyle->SetOptStat("enoum"); gPad->SetLogy()

BaseUtils::FindYield(hh, 3.5e4);

   */

  const TH1D * yy = GetHyield(hh,"findyieldtmp");
  Int_t binx = 0;
  for(Int_t ii=yy->GetMaximumBin(); ii<=yy->GetNbinsX(); ii++){
    const Double_t cont = yy->GetBinContent(ii);
    if(cont>target && target>yy->GetBinContent(ii+10)){
      binx = ii;
      break;
    }
  }
  delete yy;
  return hh->GetXaxis()->GetBinCenter(binx);
}

Double_t BaseUtils::LowMomAnchor(const Double_t xx, const Double_t anc, const Int_t opt)
{
  if(opt==0){
    const Double_t p0 = TMath::Log10(0.15);
    const Double_t p1 = TMath::Log10(5.);
    const Double_t y0 = 100;
    const Double_t y1 = anc;
    
    const Double_t aa = (y1-y0)/(p1-p0);
    const Double_t bb = y1 - aa*p1;
    
    return TMath::Min(anc, aa*xx+bb);
  }
  else if(opt==1){
    if(xx>-0.1){
      return anc;
    }
    else{
      //return TMath::Max(500. , anc/10.); //TMath::Min(anc, lpmax);
      //push it down as low as possible
      return anc/10.;
    }
  }
}

TH1D* BaseUtils::Regenerate(TH2D *hh,  Double_t anchorN, const Int_t ancOpt,   const Double_t thres)
{
  /*
  //test
gSystem->Load("../lib/libBaseUtils.so")
hh=horiginal
new TCanvas; hh->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()
hhy=BaseUtils::GetHyield(hh);

new TCanvas; hhy->Draw();  gStyle->SetOptStat("enoum"); gPad->SetLogy();

h1=(TH2D*)hh->Clone("h1")
BaseUtils::Regenerate(h1, 50, 1e2);
new TCanvas; h1->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

h1y=BaseUtils::GetHyield(h1);
new TCanvas; h1y->Draw();  gStyle->SetOptStat("enoum"); gPad->SetLogy();


h2=(TH2D*)hh->Clone("h2")
BaseUtils::Regenerate(h2, 1e4, 50);
new TCanvas; h2->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

h2y=BaseUtils::GetHyield(h2);
new TCanvas; h2y->Draw();  gStyle->SetOptStat("enoum"); gPad->SetLogy();

   */

  //if anchor too small, use the least N of all bins
  Double_t minn = 1e10;
  for(Int_t ix = 1; ix<=hh->GetNbinsX(); ix++){
    const Double_t cc = hh->Integral(ix,ix,1, hh->GetNbinsY());
    if(cc<thres){
      continue;
    }
    if(cc<minn)
      minn = cc;
  }
  if(anchorN<minn){
    anchorN = minn;
    printf("check BaseUtils::Regenerate --> anchorN modified as minimum stat %e %e\n", anchorN, minn);
  }

  //with reproducible random number with seed!=0
  TRandom3 ran(100);

  TH1D * hfityield = hh->ProjectionX("hfityield", 1, hh->GetNbinsY(),"oe");
  hfityield->Scale(0);

  for(Int_t ix=1; ix<=hh->GetNbinsX(); ix++){
    const Double_t tmpnor = hh->Integral(ix,ix,1, hh->GetNbinsY());
    const Double_t tmpx = hh->GetXaxis()->GetBinCenter(ix);
    const Double_t tmpanc = LowMomAnchor(tmpx, anchorN, ancOpt);
    if(tmpnor>tmpanc){
      //printf("test0 ix %d xx %e tmpnor %e\n", ix, hh->GetXaxis()->GetBinCenter(ix), tmpnor);

      //do not regenerate, just fake the statistics -> if not regenerated, high stat. feature not removed
      //generate tmpanc according Poisson
      Double_t ngen = 0;
      for(Int_t iy=1; iy<=hh->GetNbinsY(); iy++){
        const Double_t frac = hh->GetBinContent(ix,iy)/tmpnor;
        //http://root.cern.ch/root/htmldoc/TRandom.html#TRandom:PoissonD
        const Double_t cont = ran.Binomial((Int_t)tmpanc, frac);

        //if(cont>EPSILON){
          //printf("============= test1 iy %d yy %e n %e frac %e cont %e\n", iy, hh->GetYaxis()->GetBinCenter(iy), hh->GetBinContent(ix,iy), frac, cont);
        //}
        hh->SetBinContent(ix, iy, cont); hh->SetBinError(ix,iy,TMath::Sqrt(cont));
        ngen += cont;
      }

      hfityield->SetBinContent(ix, ngen); hfityield->SetBinError(ix, TMath::Sqrt(ngen)); 
    }
    else{
      hfityield->SetBinContent(ix, tmpnor); hfityield->SetBinError(ix, TMath::Sqrt(tmpnor));
    }
  }

  return hfityield;
}

TH1D * BaseUtils::GetHyield(const TH2D *h0, const TString hn)
{
  //only consider the dEdx in range
 //http://root.cern.ch/root/html/TH2.html#TH2:ProjectionX
  /*
    if option "e" is specified, the errors are computed.
    if option "d" is specified, the projection is drawn in the current pad.
    if option "o" original axis range of the taget axes will be
    kept, but only bins inside the selected range will be filled.
  */
  return h0->ProjectionX(hn, 1, h0->GetNbinsY(),"oe");
}

void BaseUtils::SelectXbins(TH2D *hh, const Int_t bins[])
{
  /*
//test
gSystem->Load("../lib/libBaseUtils.so")
hh=horiginal
new TCanvas; hh->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

Int_t bins[100]={170,180, 190, 320, 321,322};
h1=(TH2D*)hh->Clone("h1")
BaseUtils::SelectXbins(h1, bins);
new TCanvas; h1->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

h2=(TH2D*)hh->Clone("h2")
BaseUtils::RebinX(h2, 10, 0.6);
BaseUtils::SelectXbins(h2, bins);
new TCanvas; h2->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

   */

  //bins[ii] must be in increasing order
  Int_t ib = 0;
  for(Int_t ix=0; ix<=hh->GetNbinsX()+1; ix++){
    if(ix!=bins[ib]){
      //printf("test0 ix %d ib %d xx %e\n", ix, ib, hh->GetXaxis()->GetBinCenter(ix));
      for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
        hh->SetBinContent(ix, iy, 0); hh->SetBinError(ix,iy,0);
      }
    }
    else{
      ib++;
    }
  }
}

TH2D * BaseUtils::CutHist(TH2D *hh, const Double_t fx0, const Double_t fx1, const Double_t fy0, const Double_t fy1, const Bool_t krest/*=kFALSE*/)
{
  /*
  //test
gSystem->Load("../lib/libBaseUtils.so")
hh=horiginal
new TCanvas; hh->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

h1=(TH2D*)hh->Clone("h1")
TH2D *hr = BaseUtils::CutHist(h1, -0.5, 1.2, 20, 120, 1);
new TCanvas; h1->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()
new TCanvas; hr->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()
h1->Integral(99,99,0,10000)
h1->Integral(100,100,0,10000)
h1->Integral(442,442,0,10000)
h1->Integral(441,441,0,10000)
h1->Integral(0,100000, 40,40)
h1->Integral(0,100000, 41,41)
h1->Integral(0,100000, 242,242)
h1->Integral(0,100000, 241,241)

h2=(TH2D*)hh->Clone("h1")
hr=BaseUtils::CutHist(h2, -0.5, 1.2, 20, 120);
hr

   */

  //entries outside x0-x1, y0-y1
  TH2D *hrest = (TH2D*)hh->Clone(Form("%sout",hh->GetName()));

  //==========================================================================
  //split original, {x,y}{0,1} inclusive
  //==========================================================================
  const Int_t rx0 = hh->GetXaxis()->FindBin(fx0);
  const Int_t rx1 = hh->GetXaxis()->FindBin(fx1);
  const Int_t ry0 = hh->GetYaxis()->FindBin(fy0);
  const Int_t ry1 = hh->GetYaxis()->FindBin(fy1);
  if(rx0>=rx1 || ry0>=ry1){
    printf("BaseUtils::CutHist bad range %d %d, %d %d -- %e %e , %e %e\n", rx0, rx1, ry0, ry1, fx0, fx1, fy0, fy1);exit(1);
  }

  printf("BaseUtils::CutHist ranges [ %d %d ]*[ %d %d ]\n", rx0, rx1, ry0, ry1);

  for(Int_t ix=0; ix<=hh->GetNbinsX()+1; ix++){
    for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
      if(ix<rx0 || ix>rx1 || iy<ry0 || iy>ry1){
        hh->SetBinContent(ix,iy, 0); hh->SetBinError(ix,iy,0);
      }
      else{
        hrest->SetBinContent(ix,iy,0); hrest->SetBinError(ix,iy,0);
      }
    }
  }

  if(!krest){
    delete hrest;
    hrest=0x0;
  }

  return hrest;
}

Int_t BaseUtils::RebinX(TH2D * hh, const Int_t xrebin, const Double_t xmid, Int_t outbins[])
{
  /*
  //test
gSystem->Load("../lib/libBaseUtils.so")
hh=horiginal
new TCanvas; hh->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

h1=(TH2D*)hh->Clone("h1")
BaseUtils::RebinX(h1,10)
new TCanvas; h1->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

//total entries should be the same

h2=(TH2D*)hh->Clone("h2")
BaseUtils::RebinX(h2,10, 0.6)
new TCanvas; h2->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()
h2->GetXaxis()->GetBinWidth(319)
h2->GetXaxis()->GetBinWidth(320)
h2->GetXaxis()->GetBinWidth(321)

h3=(TH2D*)hh->Clone("h3")
BaseUtils::RebinX(h3, 7, 0.6)
new TCanvas; h3->Draw("colz"); gStyle->SetOptStat("enoum"); gPad->SetLogz()

   */

  //
  //overflow of X is not preserved in case of rebinning only for x > xmid
  //

  if(xmid <= hh->GetXaxis()->GetXmin() || xmid>=hh->GetXaxis()->GetXmax()){
    return 0;
  }
  else{
    Double_t bins[hh->GetNbinsX()+1];

    const Int_t nx = hh->GetNbinsX()+2;
    const Int_t ny = hh->GetNbinsY()+2;
    Double_t ents[nx][ny];
    for(Int_t ix=0; ix<=hh->GetNbinsX()+1; ix++){
      for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
        ents[ix][iy]=0;
      }
    }

    //============
    Int_t rmid = -999; 
    for(Int_t ix=0; ix<=hh->GetNbinsX()+1; ix++){
      if(hh->GetXaxis()->GetBinUpEdge(ix)>xmid && ix%xrebin==1){
        rmid=ix;
        break;
      }

      bins[ix]=hh->GetXaxis()->GetBinUpEdge(ix);
      for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
        const Double_t cont = hh->GetBinContent(ix,iy);
        ents[ix][iy] = cont;
      }
    }

    Int_t nbs = rmid;
    for(Int_t ix = rmid; ix<=hh->GetNbinsX(); ix+=xrebin){
      const Int_t tmp0 = ix;
      const Int_t tmp1 = ix+xrebin-1;

      bins[nbs]=hh->GetXaxis()->GetBinUpEdge(TMath::Min(hh->GetNbinsX(), tmp1));
      
      for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
        const Double_t cont = hh->Integral(tmp0, tmp1, iy, iy);
        ents[nbs][iy] = cont;
        
        //printf("filling %d %d -> %e %d %d %e\n", tmp0, tmp1, bins[nbs], nbs, iy, cont);
      }

      if(tmp1>=hh->GetNbinsX()){
        break;
      }

      nbs++;
    }

    hh->GetXaxis()->Set(nbs, bins);
    hh->Scale(0);

    Double_t totn = 0;
    Int_t outnb = 0;
    for(Int_t ix=0; ix<=nbs+1; ix++){
      if((ix<rmid && ix%xrebin==(xrebin+1)/2) || ix>=rmid){
        outbins[outnb]=ix;
        outnb++;
      }

      for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
        hh->SetBinContent(ix, iy, ents[ix][iy]);  hh->SetBinError(ix,iy,TMath::Sqrt(ents[ix][iy]));

        totn+=ents[ix][iy];
      }
    }
    hh->SetEntries(totn);

    return outnb;
  }
}

void BaseUtils::BinLogX(TAxis *axis)
{
  //
  // Method for the correct logarithmic binning of histograms
  // copied and modified from AliTPCcalibBase

  const Int_t bins = axis->GetNbins();

  const Double_t from = axis->GetXmin();
  const Double_t to = axis->GetXmax();
  if (from<EPSILON){
    //printf("BaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from);  exit(1);                                                                                                                                       
    return;
  }

  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  const Double_t factor = pow(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
    new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;
}

TH1D * BaseUtils::GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx)
{
  TH1D * h1=new TH1D(hname, "", 500, xmin, xmax);
  if(kbinlogx){
    BinLogX(h1->GetXaxis());
  }
  for(Int_t ii=1; ii<=h1->GetNbinsX(); ii++){
    const Double_t xx = h1->GetBinCenter(ii);
    const Double_t yy = func(&xx, par);
    h1->SetBinContent(ii, yy); h1->SetBinError(ii,0);
  }
  return h1;
}

TGraphAsymmErrors * BaseUtils::TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0, Char_t * sex1, Char_t * sey0, Char_t * sey1, const Bool_t kXLogToLinear, const Bool_t kprint)
{
  TTree * rawtree=(TTree*)gDirectory->Get(tname);
  if(!rawtree){
    printf("BaseUtils::TreeToGraph no tree %s\n", tname.Data());exit(1);
  }

  TTree *cuttree = rawtree->CopyTree(cut);
  TTree * tmptree=(TTree*) cuttree->Clone(Form("tmptree%s%s",tname.Data(), sy));//necessary to avoid address confusion when executed repeatedly.

  Double_t vx = -1e10, vy = -1e10, vex0 =-1e10, vex1 = -1e10, vey0 = -1e10, vey1 = -1e10;
  tmptree->SetBranchAddress(sx, &vx);
  tmptree->SetBranchAddress(sy, &vy);
  if(strlen(sex0) && strlen(sex1)){
    tmptree->SetBranchAddress(sex0, &vex0);
    tmptree->SetBranchAddress(sex1, &vex1);
  }
  if(strlen(sey0)){
    tmptree->SetBranchAddress(sey0, &vey0);
  }
  //if sey0 and sey1 are the same, vey0 will not be set from tree, has to be set by vey1
  if(strlen(sey1)){
    tmptree->SetBranchAddress(sey1, &vey1);
  }
  
  TGraphAsymmErrors * tmpgr = new TGraphAsymmErrors;
  const Int_t ntot = tmptree->GetEntries();
  for(Int_t ii=0; ii<ntot; ii++){
    tmptree->GetEntry(ii);

    //remove trivial -1e10
    if(vex0<-1e9) vex0=0;
    if(vex1<-1e9) vex1=0;
    if(vey0<-1e9) vey0=0;
    if(vey1<-1e9) vey1=0;

    //minos error low < 0 by convention
    vex0=TMath::Abs(vex0);
    vex1=TMath::Abs(vex1);
    vey0=TMath::Abs(vey0);
    vey1=TMath::Abs(vey1);

    if(strcmp(sey0,sey1)==0){
      vey0=vey1;
    }

    if(kXLogToLinear){
      const Double_t xtrue = pow(10, vx);
      const Double_t xlow = pow(10, vx-vex0);
      const Double_t xhigh = pow(10, vx+vex1);
      vx = xtrue;
      vex0 = xtrue - xlow;
      vex1 = xhigh - xtrue;
    }

    tmpgr->SetPoint(ii, vx, vy);
    tmpgr->SetPointError(ii, vex0, vex1, vey0, vey1);
    if(kprint){
      printf("BaseUtils::TreeToGraph %s %s %s -- %d %e %e -- %e %e %e %e\n", tname.Data(), sx, sy, ii, vx, vy,  vex0, vex1, vey0, vey1);
    }
  }
  delete tmptree;

  TGraphAsymmErrors * gr = (TGraphAsymmErrors *)tmpgr->Clone(Form("gr%s%s",tname.Data(), sy));
  delete tmpgr;
  return gr;
}

TH1D* BaseUtils::ToPDF(const TH1 *hraw, const TString hn)
{
  TH1D * hist = (TH1D*) hraw->Clone((hn+hraw->GetName())+"pdf");
  Double_t ntot = 0;
  hist->Scale(0);

  //only consider region in range, ignore over/underflow. Important for fitting experiment data where under/overflow is ignored.
  const Int_t x0 = hist->GetXaxis()->GetFirst();
  const Int_t x1 = hist->GetXaxis()->GetLast();

  for(Int_t ii=x0; ii<=x1; ii++){
    const Double_t cont = hraw->GetBinContent(ii);
    const Double_t err = hraw->GetBinError(ii);
    //skip empty bins
    if(cont<EPSILON){
      continue;
    }

    //check scale
    //EPSILON to strong
    const Double_t eps = 1e-3;
    if(fabs(TMath::Sqrt(cont)-err)>eps){
      printf("BaseUtils::ToPDF is scaled!! %e %e %e\n", cont, err, fabs(TMath::Sqrt(cont)-err)); exit(1);
    }

    const Double_t nn = cont;
    hist->SetBinContent(ii,nn); hist->SetBinError(ii, TMath::Sqrt(nn));
    ntot += nn;
  }

  if(ntot<EPSILON){
    printf("BaseUtils::ToPDF ntot<epsilon ! %f %f\n", hist->GetEntries(), ntot);
    exit(1);
  }
  
  const Double_t bw = hist->GetBinWidth(1);
  for(Int_t ib=x0; ib<=x1; ib++){
    const Double_t cont = hist->GetBinContent(ib);
    if(cont<EPSILON)
      continue;

    const Double_t eff = cont/ntot;
    const Double_t pdf = eff/bw;

    //const Double_t dpdf = sqrt(eff*(1-eff)/ntot) / bw;
    //use Poisson error! Relevant for fitting the PDF with binned likelihood
    const Double_t dpdf = TMath::Sqrt(cont)/ntot/bw;
    hist->SetBinContent(ib, pdf); hist->SetBinError(ib, dpdf);
  }

  hist->SetEntries(ntot);

  return hist;
}


TH2D* BaseUtils::NormalHist(const TH2D *hraw, const Double_t thres, const Bool_t kmax)
{
  //
  //normalized histogram
  //

  TH2D *hh=(TH2D*)hraw->Clone(Form("%snor",hraw->GetName()));
  hh->Scale(0);

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  Double_t hmax = -1e10;
  Double_t hmin = 1e10;
  Double_t nent = 0;
  for(Int_t ix=x0; ix<=x1; ix++){

    //if option "e" is specified, the errors are computed. if option "o" original axis range of the taget axes will be kept, but only bins inside the selected range will be filled.
    TH1D * sliceh = hraw->ProjectionY(Form("tmpnormalhist%sx%d", hh->GetName(), ix), ix, ix, "oe");
    const Double_t tot = sliceh->Integral(1, sliceh->GetNbinsX());//bad: sliceh->GetEntries();

    //printf("test0 ix %d xx %e tot %e int %e\n", ix, hraw->GetXaxis()->GetBinCenter(ix), tot, sliceh->Integral(0,100000));

    TH1D * pdfh=0x0;

    if(tot>EPSILON){
      nent += tot;
    
      Double_t imax = -999;

      if(!kmax){
        pdfh = ToPDF(sliceh,"tmp");
      }
      else{
        imax = sliceh->GetBinContent(sliceh->GetMaximumBin());
      }

      for(Int_t iy=y0; iy<=y1; iy++){
        const Double_t cont = kmax ? sliceh->GetBinContent(iy)/imax : pdfh->GetBinContent(iy);
        const Double_t ierr = kmax ? sliceh->GetBinError(iy)/imax   : pdfh->GetBinError(iy);
        if(tot>thres && cont>0){
          //printf("test0 ix %d xx %e iy %d yy %e tot %e cont %e\n", ix, hraw->GetXaxis()->GetBinCenter(ix), iy, hraw->GetYaxis()->GetBinCenter(iy), tot, cont);
          hh->SetBinContent(ix, iy, cont); hh->SetBinError(ix,iy, ierr);
          if(cont>hmax) hmax = cont;
          if(cont<hmin) hmin = cont;
        }
      }
    }

    delete pdfh;
    delete sliceh;
  }

  hh->SetEntries(nent);
  hh->SetMinimum(0.99*hmin);
  hh->SetMaximum(1.1*hmax);
  return hh;
}

Double_t * BaseUtils::GetAxisArray(TAxis * aa)
{
  TArrayD * xs = new TArrayD(*(aa->GetXbins()));
  if(xs->GetSize()==0){
    const Int_t nbin=aa->GetNbins();
    xs->Set(nbin+1);
    for(Int_t ii=0; ii<=nbin; ii++){
      xs->SetAt(aa->GetBinUpEdge(ii),ii);
    }
  }

  Double_t * bins = xs->GetArray();
  if(!bins){
    printf("BaseUtils::GetAxisArray bins null!! %d\n", xs->GetSize());exit(1);
  }
  return bins;
}
